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ABSTRACT 

NASHUA is a coupled finite element/boundary element capability built around 
NASTRAN for calculating the low frequency far-field acoustic pressure field radiated or 
scattered by an arbitrary, submerged, three-dimensional, elastic structure subjected to either 
internal time-harmonic mechanical loads or external time-harmonic incident loadings. This 
paper describes the formulation and use of NASHUA for solving such structural acoustics 
problems when the structure is fluid-filled. NASTRAN is used to generate the structural 
finite element model and to perform most of the required matrix operations. Both fluid 
domains are modeled using the boundary element capabililty in NASHUA, whose matrix 
formulation (and the associated NAS I RAN DMAP) for evacuated structures can be used 
with suitable interpretation of the matrix definitions. Alter computing surface pressures and 
normal velocities, far-field pressures are evaluated using an asymptotic form of the 
Helmholtz exterior integral equation. The proposed numerical approach is validated by 
comparing the acoustic field scattered from a submerged fluid-filled spherical thin shell to 
that obtained with a series solution, which is also derived in this paper. 


INTRODUCTION 

Two basic problems in computational structural acoustics are (1) the calculation of the 
acoustic pressure field radiated by a general submerged three-dimensional elastic structure 
subjected to internal time-harmonic loads, and (2) the calculation of the acoustic pressure 
scattered by such a structure subjected to an incident time-harmonic wavetrain. The most 
common, as well as the most accurate, approach for solving these problems at low 
frequencies is to couple a finite element model of the structure with a boundary element 
model of the surrounding fluid. This is the approach taken by NASHUA, which is a 
boundary element program built around NASTRAN, a widely-used finite element computer 
program for structural dynamics. 

Several previous papers (Ref. 1-4) described the basic formulation and development for 
acoustic radiation and scattering from evacuated structures. Here we describe the 
formulation and use of NASHUA for modeling submerged structures which are fluid-filled. 
Internal fluid can occur because the structure is Iree-llooded or contains fluid-filled tanks. It 
is possible to use existing NASTRAN capability to model the interior fluid with finite 
elements (Ref. 5-7), but three-dimensional models with large numbers of fluid degrees of 
freedom might result. An attractive alternative to the fluid finite element model is to 
represent the contained fluid using a boundary element approach. In principle, any computer 
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program capable of generating the appropriate boundary element matrices for an exterior 
fluid is also capable of generating such matrices for the complementary region (the interior 
region). NASTRAN’s versatility in user-controlled matrix operations (DMAP) makes the 
implementation of such an approach straightforward. 


THEORETICAL APPROACH 

The basic theoretical development for NASHUA’S radiation and scattering approach 
for evacuated structures has been presented in detail previously (Ref. 1-4). For 
completeness, this paper summarizes that approach and describes the changes necessary to 
model the interior fluid with boundary elements in the same procedure. There is no 
requirement that the interior and exterior fluids be the same. 


The Surface Solution for Evacuated Structures 

Consider any submerged three-dimensional, evacuated elastic structure subjected to 
either internal time-harmonic loads or an external time-harmonic incident wavetrain. If the 
structure is modeled with finite elements using NASTRAN, the resulting matrix equation of 
motion can be written as 

Zv = F-GAp, 0 ) 

where matrix Z (of dimension s x s) is the structural impedance, vector v (s x r) is the 
complex velocity amplitude for all structural DOF (wet and dry) using the coordinate systems 
selected by the user, vector F (s x r) is the complex amplitude of the mechanical forces 
applied to the structure, matrix G (s x f) is the rectangular transformation of direction 
cosines to transform a vector of outward normal forces at the wet points to a vector of 
forces at all points in the coordinate systems selected by the user, matrix A (f x f) is the 
diagonal area matrix for the wet surface, and vector p (f x r) is the complex amplitude of 
total pressures (incident + scattered) applied at the wet grid points. In this equation, the 
time dependence exp(iwt) has been suppressed. In the above dimensions, s denotes the total 
number of independent structural DOF (wet and dry), f denotes the number of fluid DOF 
(wet points), and r denotes the number of load cases. In general, the surface areas, the 
normals, and the transformation matrix G are obtained in NASHUA from the NAS IRAN 
calculation of the load vector resulting trom an outwardly directly static unit pressure load on 
the structure’s wet surface. 

In Eq. 1, the structural impedance matrix Z, which converts velocity to force, is given 
by 

Z = (-arM +iwB + K)/(iw), ( 2 ) 

where M, B, and K are the structural mass, viscous damping, and stiffness matrices, 
respectively, and u is the circular frequency of excitation. For structures with a nonzero loss 
factor, K is complex. In addition, K can include the differential stiffness effects of 
hydrostatic pressure, if any (Ret. 3). A standard NASTRAN finite element model of the 
structure supplies the matrices K, M, and B. 

For the exterior fluid domain, the total fluid pressure p satisfies the Helmholtz 
differential equation 
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+ k 2 p = o, 


( 3 ) 


where k = cc/c is the acoustic wave number, and c is the fluid sound speed. Equivalently, for 
the exterior fluid, p is the solution of the Helmholtz integral equations (Ref. 8) 


/ s P( x ) dS - / s q(x)D(r)dS = 


p(x')/2 — pj, x ' on 

P( x ')-Pi> x' in E, 

-n,. x' in I, 


( 4 ) 


where S, E, and I denote the surface, exterior, and interior domains, respectively, pj is the 
incident free-field pressure (if any), r is the distance from x to x' (Fig. 1), D is the frce-space 
Green’s function 
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Fig. 1. Notation for Helmholtz Integral Equation 


p is the fluid mass density, and v n is the outward normal velocity on S. As shown in Fig. 1, x 
in Eq. 4 is the position vector for a typical point Pj on the surface S, x' is the position vector 
for the point P; on the surface or in the exterior field, the vector r = x' - x, and n is the unit 
outward normal at Pj. We denote the lengths of the vectors x, x', and r by x, x', and r, 
respectively. The normal derivative of the Green’s function D is (Ref. 1) 


dD(r) 

On 


—ikr 


— (ik + — ) cos /?, 

4 Trr r 


( 7 ) 


where 0 is the angle between the normal n and the vector r, as shown in Fig. 1. 

All three integral equations in Eq. 4 are needed for exterior fluids. The surface 
equation provides the fluid impedance at the fluid-structure interface. Since the surface 
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equation exhibits non-uniqueness at certain discrete characteristic frequencies (Ref. 9), the 
interior equation is used to provide additional constraint equations which ensure the required 
uniqueness. The exterior equation is used to compute the exterior pressure field once the 
surface solution (which includes the fluid pressure and its gradient) is known. 


The substitution of Eqs. 6 and 7 into the surface equation (4) yields 

y jv ikr 1 

J2i2Ll _ jf p (x) (ik + y) cos p dS = ice/> £ v„(x) -^-dS + Pi’ x ' on S - ( 8) 


This integral equation relates the total pressure p and normal velocity v n on S. II the 
integrals in Eq. 8 are discretized for numerical computation (Ref. 1), we obtain the matrix 
equation (for the exterior fluid) 

Ep = Cv n + p, , W 


where vector p (ol dimension f x r) is the vector of complex amplitudes ol the total pressure 
on the Structure’s wet surface, matrices E and C (both f x f) are fully-populated, complex, 


on the structure’s wet surface, matrices E and C (both 1 x f) are lully-populated, ec 
nonsymmetric, and frequency-dependent, and vector pj (f x r) is the complex amplitude 
the incident pressure vector. The number of unknowns in this system is f, the number 
wet points on the lliiid-structurc interface. 


of 


The normal velocities v„ in Eq. 9 are related to the total velocities v by the same 
rectangular transformation matrix G: 

v„ = G T v, 00) 

where T denotes the matrix transpose. If velocities v and v n are eliminated from Eqs. 1, 9, 
and 10, the resulting equation for the coupled fluid-structure system is 


(E + CG t Z _ 1 GA) p = CG T Z _1 F + pj. (ID 

This equation is solved for the total surface pressures p, since the rest of the equation 
depends only on the geometry, the material properties, and the frequency. Since the two 
right-hand side terms in Eq. 1.1 correspond to mechanical and incident loadings, only one ol 
the two terms would ordinarily be present lor a given case. I he details ol the incident 
pressure vector p! for scattering problems were presented previously (Ref. 2). and will not be 
repeated here. 


The velocity vector v for all structural DOF is recovered by solving Eq. 1 for v: 


v = Z -1 F — Z _1 GAp. ( 12 > 

The surface normal velocity vector v n is recovered by substituting this solution for v into Eq. 

10 . 


Modeling Interior Fluid 

The theoretical development presented in the preceding section can be modified slightly 
to account also for an interior fluid. The wave equation, Eq. 3, applies also to interior fluids. 
Although all three integral equations in Eq. 4 are generally needed for exterior fluids, only 
the surface equation is needed to represent the surface impedance of interior fluids. Eq. 4a 
also applies to interior fluids if the incident pressure p! is set to zero, and the normal vector 
n is negated. That is, the surface integral equation applies to both exterior and interior fluids 
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so long as the unit normal is always directed from the structure into the fluid. One other 
consideration, perhaps unique to NASHUA, is that wet surface curvatures (which are 
needed in the calculation of the "self" terms in matrix E) are negative at some interior points 
(Ref. 1). 

A matrix equation similar to Eq. (9) is therefore obtained for the interior fluid except 
that the incident pressure p t is zero. The fluid matrices E and C are different for exterior 
and interior domains (even if the separating surface S has infinitesimal thickness) because 
the normals are of opposite sign. 


Two-Fluid Formulation 

Denote the exterior fluid as Fluid 1 and the interior fluid as Fluid 2, and use the 
subscripts 1 and 2 to refer to these two domains. Also define new pressure and normal 
velocity unknowns p and v n which include the solutions for both fluid domains: 



Since there is no direct fluid coupling between the interior and exterior fluids, and the 
incident pressure vanishes in the interior domain, Eqs. 1, 9, 10, and 11 apply also to the 
two-fluid problem if the new definitions in Eq. 13 are used, and the matrices A, G, E, C, 
and pi are re-defined as 

p. = { P 4 oo 

The principle benefit of formulating the two-fluid problem in this way is that the required 
modifications to extend the procedure to three or more independent fluid domains is now 
clear. 



The Far-Field Calculation 

With the solution for the total pressures and velocities on the surface, the exterior 
Helmholtz integral equation, Eq. 4b, can be integrated to obtain the radiated (or scattered) 
pressure at any desired location x' in the exterior field. We first substitute Eqs. 5-7 into Eq. 
4b to obtain 

1 e” ikr 

p(x') = / [iwy>v n (x) + (ik + — )p(x) cos f3\ — dS, x' in E. (IS) 

s r 47rr 

In applications, however, the field pressures generally of interest are in the far-field, so we 
use instead the asymptotic form of Eq. 15 (Ref. 1): 

lkp ikx 

p(x') = 47rx , J s [pcv n (x) -F P(x) cos ^]e ,kxcosa dS, x' in E, x' » d, (16) 

where d is a characteristic dimension of the structure, and a is the angle between the vectors 
x and x' (Fig. 1). For far-field points, cos /? is computed using the asymptotic approximation 


154 



a X 

cos p — ► n*— 
x 


(17) 


For both Eqs. 15 and 16, numerical quadrature is used. 


OVERVIEW OF SOLUTION PROCEDURE 

The NASHUA solution procedure uses NASTRAN to generate the matrices K, M, B, 
and F and to generate sufficient geometry information so that the matrices E, C, G, A, and 
Pi can be computed by a separate program called SURF. Then, NASTRAN DMAP is used 
to form the matrices appearing in Eq. 11, which is solved for the total pressures p (in both 
fluid domains) using the block solver OCSOLV (Ref. 10). Next, NASTRAN DMAP is used 
to recover the surface normal velocities v n and the vector v of velocities at all structural 
DOF (NASTRAN’s "g-set"). This step completes the surface solution. Then, with the total 
pressures and velocities on the (exterior) surface, the asymptotic (far-ficld) form of the 
Helmholtz exterior integral equation is integrated in program FAROUT to compute the far- 
field radiated pressures. Various tables and graphical displays are generated. 

The overall setup of the solution procedure is organized into four steps. In Step 1, a 
separate NASTRAN structural model is prepared and run for each unique set of symmetry 
constraints and each fluid region. Since, for general three-dimensional analysis, up to three 
planes of reflective symmetry are allowed, there would be one, two, four, or eight such runs 
for each fluid region. Since the purpose of this step is to generate a file containing geometry 
information and a checkpoint file for subsequent use in the other steps, the only difference 
between the two runs associated with a given symmetry case is the specification of the 
outwardly directed unit pressure load which defines the wet surface for a given fluid region. 

For each symmetry case and drive frequency, several programs are run sequentially to 
form Step 2. For each fluid region, the SURF program reads the geometry file generated by 
NASTRAN in Step 1 and, using the Helmholtz surface and interior integral equations, 
generates the fluid matrices E 1? E 2 , Cj, and Co, the area matrices Aj and A 2 , the structure- 
fluid transformation matrices Gj and G 2 , the incident pressure vector pn, and a geometry 
file to be used later by the far-field integration program FAROUT in Step 3. In addition, a 
partitioning vector is generated to facilitate the merging and partitioning of the various 
matrices associated with the two fluid domains. 

The two SURF jobs in Step 2 are followed by a NASTRAN job which takes the 
structural matrices K, M, B, and F from Step 1 and the matrices generated by the SURF 
jobs and forms the matrices in Eq. 11, where the definitions in Eq. 14 apply. Eq. 11 is then 
solved for the total surface pressure vector p by program OCSOLV, which is a general out- 
of-core block solver designed specifically for large, full, complex, nonsymrnetric systems of 
linear, algebraic equations. NASTRAN is then re-entered in Step 2 with p so that the 
velocities v and v n can be recovered using DMAP operations. The surface pressures, 
normal velocities, and full g-set displacements are then reformatted, sorted, and merged into 
a single file (for each symmetry case) using program MERGE. Recall that there are one, 
two, four, or eight possible symmetry cases. 

Steps 1 and 2 are repeated for each symmetry case. After all symmetry cases have 
been completed and merged, program FAROUT (Step 3) combines the symmetry cases and 
integrates over the surface. The far-field pressure solution is obtained by integrating the 
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surface pressures and velocities using the asymptotic (far-held) form of the exterior 
Helmholtz integral equation, Eq. 16. Output from FAROUT consists of both tables and files 
suitable for various types of plotting. 

The remaining steps in the NASHUA procedure are for graphical display. Deformed 
structural plots of the frequency response are obtained by restarting NASTRAN (Step 4) 
with the checkpoint hie from Step 1 and a results file from FAROUT. In addition, animated 
plots can be generated on the Evans & Sutherland PS-330 graphics terminal using the 
CANDI program written for the DEC/VAX computer by R.R. Lipman of DTRC (Ref. 11). 
X-Y plots of various quantities (both surface and far-held) versus frequency may be obtained 
using IPLOT or other interactive plotting programs (Ref. 12). Polar plots of the far-held 
sound pressure levels in each of the three principal coordinate planes can also be generated 
using the interactive graphics program FAFPLOT (Ref. 1). 


DMAP ALTERS 

Several DMAP alters arc used in the overall NASHUA procedure to implement the 
procedure described in preceding section. For Step 1, the following alter is used: 


ALTER 

ALTER 

ALTER 

GP3 

ALTER 

SSG'l 


SSG2 

OUTPUT2 

OUTPUT2 

OUTPUT2 

PAR A MR 

COND 

PARAMR 

DIAGONAL 

ADD 

RBMG2 

SSG3 

SDR1 

TA1 

DSMG1 

EQUIV 

COND 

MCE2 

LABEL 

EQUIV 

COND 


1 S NASHUA STEP 1, COSMIC 1988 RF8 (REVISED 12/7/89) 

2,2 $ DELETE PRECHK 
21,21 S REPLACE GP3 

GEOM3,EQEX1N,GEOM2/SLT,GPTT/S,N,NOGRAV/NEVER=1 $ 
117,117 $ REPLACE FRRD 

SLT.BGPD r,CSTM,SIL,EST,MPT,GPTT,EDT,MGG,CASECC,DIT,/ 
PG „ , ,/LUSET/NSKIP $ PG 

USET,GM,YS,KFS,GO,DM,PG/QR,PO,PS,PL $ PL 
A XI C , BG PD 1, EQ EX I N , U SET, PG $ 

PL , C STM , ECT , , $ 

„„//- 9 $ 

//*EQ*//C,Y,IISP=0./0.////NOHSP $ 

LBL4D,NOIISP $ SKIP DIFF. STIFF. IF NO HYDRO. P 
//*COMPLEX//C,Y,IISP=0./0./HSPC $ IISP+PO 
KAA/KDIAG/*SQUARE*/1.0 $ 

KAA,KDIAG/KAAD/(1 .0,0.0)/(l.E-6,0.) $ 

KAAD/LLL $ FACTOR KAA 

LLL,KAAD,PL,LOO,KOO,PO/ULV,UOOV,RULV,RUOV/OMl'I7 
V,Y,IRES=-1/1/S,N,EPSI $ STATIC SOLUTION 
USET,PG,ULV,UOOV,YS,GO,GM,PS,KFS,KSS,/UGV,PGG,QG/l/ 
*BKL0* $ RECOVER DEPENDENT DISPLACEMENTS 
ECT,EPT,BGPDT,SIL,GPTT,CSTM/X1,X2,X3,ECPT,GPCT/LUSET/ 
NOSIMP/O/NOGENL/GENEL $ TABLES FOR DIFF. STIFF. 
CASECC,GPTr,SIL,EDT,UGV,CSTM,MPT,ECPT,GPCT,DlT/KDGG/ 
S,N,DSCOSET $ DIFF. STIFF. MATRIX 

KDGG,KDNN/MPCF2 / MGG.MNN/MPCF2 $ EQUIV IF NO MFC’S 
LBL1D.MPCF2 $ TRANSFER IF NO MFC’S 
USET,GM,KDGG,„/KDNN,„ $ MFC’S ON DIFF. STIFF. 

LBL1D $ 

KDN N , K DFF/S I N GLE / MNN.MFF/SINGLE S EQUIV. IF NO SPC’S 
LBL2D, SINGLE S TRANSFER IF NO SPC’S 
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SCE1 

LABEL 

EQUIV 

COND 

SMP2 

LABEL 

PAR A MR 

ADD 

ADD 

EQUIV 

LABEL 

DIAGONAL 

ADD 

ADD 

FRRD 


CHKPNT 
CHKPNT 
EXIT $ 

END ALTER $ 


USET,KDNN,,,/KDFF,KDFS,KDSS,,, $ SPC’S AND DIFF. STIFF. 
LBL2D $ 

KDFF.KDAA/OMIT / MFF,MAA/OMIT $ EQUIV. IF NO OMITS 
LBL3D,OMIT $ TRANSFER IF NO OMITS 
USET,GO,KDFF/KDAA $ OMITS AND DIFF. STIFF. 

LBL3D $ 

//*SUBC*////MHSPC//HSPC $ NEGATE HYDRO. P 
KDD,KDAA/NEWKDD/(1.0,0.0)/MHSPC $ 
KFS,KDFS/NEWKFS/(1.0,0.0)/MFISPC $ 
NEWKDD,KDD//NEWKFS,KFS $ 

LBL4D $ END OF DIFF. STIFF. EFFECTS (HSP) 
KDD/IDENT/*SQUARE*/0. $ D-SET IDENTITY 
IDENT,/IDM/(1. 0,0.0) $ ANOTFIER D-SET IDENTITY 
IDENT,/ZERO/ (0.0, 0.0) S D-SET ZERO MATRIX 
CASEXX,USETD,DLT,FRL,GMD,GOD,IDENT,ZERO,IDM,,DIT/ 
UDVF,PSF,PDF,PPF/*DISP*/*DIRECT*/LUSETD/MPCF1/ 
SINGLE/OMIT/NONCUP/FRQSET $ PDF, KDD=I, BDD=0, MDD=i 
MDD,KDD,BDD,PDF,PSF,PPF,EQDYN,USETD,GOD,GMD $ 
KFS,BGPDT,ECT,EQEXIN,GPECT,SIL $ 


The above alter does not depend on whether the fluid is interior or exterior to the structure. 
The Step 2 alters, however, depend on whether an interior fluid is present. For Step 2A, the 
following alter is used: 


ALTER 1 S NASHUA STEP 2A, COSMIC 1988 RF8 (REVISED 11/7/89) 

ALTER 2,167 $ REPLACE ALL AFTER ’BEGIN’ AND BEFORE ’END’ 

INPUTT2 /DAT2,„,//13 $ INTERNAL FLUID 

INPUTT2 /DAT,,„//11 $ READ SURF MATRIX FROM UT1 

MATPRN DAT,DAT2,„ $ 

PARAML DAT//*DMI*/l/8/RIGD $ GET RIGID FLAG 

PARAMR //*EQ*//RIGD/0.////ELAST $ SET ELAST=-1 IF ELASTIC 
COND LBL9D,ELAST $ IF ELASTIC, JUMP OVER RIGID/SOFT 

PARAMR //*EQ*//RIGD/2.////SOFT $ SET SOFT=-l IF SOFT BD. 

COND LBL9E.SOFT $ IF SOFT BOUNDARY, JUMP OVER RIGID 

INPUTT2 /E,PI,VEKC„//11 $ READ SURF MATRICES FROM UT1 

OUTPUT2 PI,E,„ //-l $ INPUTT2 FILE IS OVER WRITTEN (UT1) 
OUTPUT2 „„ //- 9 $ EOF 

CHKPNT DAT,VEKC $ 

EXIT $ 

LABEL LBL9E $ BEGINNING OF SOFT ANALYSIS 

INPUTT2 /CT,PI,VEKC„//11 $ READ SURF MATRICES FROM UT1 

TRNSP CT/C $ 

ADD PI, /MPI/(-l .0,0.0) $ NEGATE PI 

OUTPUT2 MPI,C,,, //-l $ INPUTr2 FILE IS OVER WRITTEN (UT1) 
OUTPUT2 „„ //-9 $ EOF 

CHKPNT DAT,VEKC $ 

EXIT $ 

LABEL LBL9D $ BEGINNING OF ELASTIC ANALYSIS 

INPUTT5 /G2,A2,„//14 $ INTERNAL FLUID 
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INPUTT2 

INPUTT5 

INPUTT2 

MATPRN 

MERGE 

MERGE 

MERGE 

MERGE 

MERGE 

MATPRN 

PARAME 

PARAMR 

PARAMR 

PARAMR 

PARAMR 

PARAMR 

ADD5 

MPYAD 

DECOMP 

IBS 

IBS 

ADD 

ADD 

MPYAD 

MPYAD 

MPYAD 

MPYAD 

MERGE 

MERGE 

MERGE 

EQUIV 

MERGE 

MERGE 

MERGE 

ADD 

EQUIV 

OUTPUT2 

OUTPUT2 

CIIKPNT 

CIIKPNT 

END ALTER 


/C2T,E2,PI2,VEKC2,//13 $ INTERNAL FLUID 
/G1,A1,,,//12 $ READ SURF MATRICES FROM UT2 
/C1T,E1,PI1,VEKC,FVEC//11 $ READ SURF MATRICES 
F VEC , , , , $ 

A1„,A2,FVEC,/A/-1 $ 

El, , ,E2,FVEC,/E/-1 $ 

CIT,,,C2T,FVEC,/CT/-1 $ 

Gl„G2„FVEC,/G/0 S 
PIl„„,FVEC/PI/0 $ 

VECM.VECS,,, $ 

DAT//*DMI*/l/2/FREQ $ GET FREQ FROM DAT 
//*COMPLEX*//FREQ/0./FREQC $ FREQ+DO 
//*MPYC*////W/FREQC/(6.283185,0.) $ OMEGA 
//*MPYC*////IW/W/(0.,1.) $ DOM EGA 
//*MPYC*////WW/W/W $ OMEGA**2 
//*SUBC*////MWW//WW $ -OMEGA**2 
MDD,KDD,BDD,,/Y/MWW/(1.0,0.0)/1W $ 

G.A./GA/O $ 

Y/L,U/1//S,N,MINDIAG///S,N,SING $ 

L,U,GA/YIGA/1 S 
L,U .PDF/YIF/ 1 $ 

YIGA./ZIGA/1W S 
YIF./ZIF/IW $ 

G,ZIGA,/GTZIGA/1 $ 

CT,GTZ1GA,E/I1/1 $ LHS 
G,ZIF,/GTZIF/1 $ 

CT,GTZIF,/Q/1 S MECHANICAL RHS 
DUM,,PDF,,VECM,/PDF1/1 $ MERGE IN 0 COLUMNS 
DUM,,PSF,,VECM,/PSI‘T/1 $ MERGE IN 0 COLUMNS 
DUM,,PPF,,VECM,/PPF1/1 $ MERGE IN 0 COLUMNS 
PDF - LPDF//PSF1 ,PSF//PPF1,PPF $ 

DUM,,Q,,VECM,/RIIS1/1 $ MERGE IN ZERO COLUMNS 
DUM,,GTZIF,,VECM,/G rZIFE/1 $ MERGE IN 0 COLUMNS 
DUM,,PI,,VECS,/RTIS2/1 $ MERGE IN ZERO COLUMNS 
RH S 1 , RHS2/RHS $ ADD MECH. AND INC. RHS 
USETD,DUMl//GOD,DUM2//GMD,DUM3//KFS,DUM4 $ 
RIIS,II„, HA $ INPU'm FILE IS OVER WRITTEN (UT1) 

II -9 $ EOF 

GTZIGA,GTZIFE,GA,PDF,L,U,PSE,DAT,VEKC,FVEC $ 
USE I D,GOD,GMD,KES $ 

$ 


The differences between this alter and one used for submerged evacuated structures are due 
to the need to read and combine two sets of SURF matrices, one for each fluid domain. For 
Step 2B, the following alter is used: 


ALTER 1 $ NAvSHUA STEP 2B, COSMIC 1988 RF8 (REVISED 11/7/89) 

ALTER 2,167 $ REPLACE ALL AFTER ’BEGIN’ AND BEFORE ’END’ 

INPUTT2 /PC,,„//11 $ READ PRESSURES FROM BLOCK SOLVER (UT1) 
PAR I N PC„FVEC/Pl„,/0 S REMOVE INTERNAL FLUID DOF 

PARTN Pl,,VEKC/P,„/0 $ REMOVE CHIEF DOF FROM P 
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COND 

OUTPUT2 

OUTPUT2 

MATPRN 

EXIT 

LABEL 

MPYAD 

MPYAD 

FBS 

SDR1 

PARTN 

PARTN 

OUTPUT2 

OUTPUT2 

MATPRN 

ENDALTER 


LBL9D, BLAST $ IF ELASTIC, JUMP OVER RIGID/SOFT 
DAT,P,„ //-l $ INPUTT2 FILE IS OVER WRITTEN (UT1) 

„„ //-9 $ EOF 

DAT,P,„ $ FOR SOFT BOUNDARY, P REPRESENTS VELOCITY 
$ 

LBL9D $ ELASTIC ANALYSIS 

GTZIG A , PC , GTZIFE/ VNC/0/-1 S NORMAL VELOCITIES 
GA,PC,PDF/FA/0/-l S A-SET FORCES 
L,U,FA/UDVF/1 $ A-SET DISPLACEMENTS 
USETD,,UDVF„,GOD,GMD,PSF,KFS,,/UPVC,,QPC/l/ 
^DYNAMICS* $ 

VNC„FVEC/Vl,,,/0 $ REMOVE INTERNAL FLUID DOF 
Vl„VEKC/VN,,,/0 $ REMOVE CHIEF DOF FROM VN 
DAT,P,VN,UPVC, //-l $ INPUTT2 FILE IS OVER-WRI ITEN 
,,„ H-9 $ EOF 
DAT,P,VN,, $ 

$ 


This alter differs from one for evacuated structures because of the presence of several matrix 
partitionings to remove the internal fluid DOF from the solution vectors before the solutions 
are merged with the results for other frequencies. 


NUMERICAL EXAMPLE 

Here we illustrate and validate the two-fluid boundary element formulation developed 
above by solving the problem of acoustic scattering from a submerged fluid-filled spherical 
thin shell. The incident loading is a time-harmonic planar wavctrain, as shown in Fig. 2. 
The specific problem solved has the following characteristics: 


shell mean radius (a) 

5 m 

shell thickness (h) 

0.15 m 

shell Young’s modulus (E) 

2.07 x 10 u N/m 2 

shell Poisson’s ratio (//) 

0.3 

shell density ( p s ) 

7669 kg/nF 

shell loss factor (//) 

0.01 

lluid density ( p ) 

1000 kg/m 3 

fluid sound speed (c) 

1524 m/s 


The same lluid is used for both the exterior and interior fluid domains. The solution of this 
problem exhibits rotational symmetry about the spherical axis parallel to the direction of 
wave propagation. The benchmark solution to which the numerical results will be compared 
is a series solution, the derivation of which is summarized in the next section. 


Series Solution 

The series solution for scattering from a submerged evacuated spherical thin shell was 
presented by Junger and Feit (Ref. 13). Here we summarize that solution and indicate the 
modification necessary to include the addition of an interior lluid which fills the spherical 
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Fig. 2. Plane Wave Scattering from a Fluid-Filled Spherical Shell 

volume. 


In general, the series solution for plane wave scattering from a submerged, evacuated, 
spherical thin shell involves computing the impedances of the shell and exterior fluid, the 
scattered field due to rigid body effects, and the radiated field due to elastic (vibrational) 
effects. The shell impedance (the ratio of pressure to normal velocity) for the nth 
axisymmetric shell mode is 


ipsc p h [^-(o^) 2 ] [rf-(n^) 2 ] 

^ a [fl” — (1+/?“)(V+X n — 1] 


(18) 


where p s is the structural mass density, c p = Ve/[/) s ( 1-^)] , E is Young’s modulus, v is 
Poisson’s ratio, 0= tea/c p is tha_ dimensionless frequency, h is the shell thickness, a is the 
shell mean radius, ^ = h/(avl2), and X n = n(n-HL). The quantities and are the 
upper and lower shell resonance dimensionless frequencies, respectively. They are the 
solutions of the characteristic equation 


if - [ 1 +3 //+X n —{? ( 1 — //— X n — tA n ] if 

+ (\i — 2)(l— i/*')+/?‘'[X n J ’-4X n “+\(5— if)— 2(1— /r - )] = 0. (19) 


The impedance of the exterior fluid, found by using the Green’s function and identity for the 
exterior fluid, is 


z 


n 


\pc 


h n (ka) 
h' n (ka) ’ 


( 20 ) 


where h n is the Bessel’s function of the third kind of order n. 


Thus, Junger and Feit showed that the far-field scattered pressure is 
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pOM)= 


, R>>a, 


( 21 ) 


ie ikR Po (2n+l)P n (cosg) pc 

kR ,h h 'n(ka) [ j ,A aj (Z n +z n )(ka) 2 h' n (ka) 

where R is the distance to the field point, 0 is the angle from the z-axis, p G is the incident 
pressure, P n is the Legendre polynomial of order n, and j n is the Bessel’s function of the first 
kind of order n. The two terms in the bracketed expression correspond to rigid body and 
radiated effects, respectively. 


The above expression for the pressure scattered from an evacuated shell can be 
extended to include the effects of the interior fluid merely by replacing the exterior fluid 
impedance z n in Eq. 21 with the sum of the fluid impedances for the exterior and interior 
fluids. It can be shown, by using the Green’s function and identity for the interior domain, 
that the interior impedance, denoted f n , is given by 


Sn 


— ipc 


jn( ka ) 

j'n(ka) ‘ 


( 22 ) 


We note the resemblance between Eqs. 20 and 22 for the exterior and interior domains, 
respectively. 

The computer program used to evaluate this series solution is a modification of a 
program called SCATSPHERE written by F.M. Henderson, a retired employee of DTRC. 
SCATSPHERE in turn is a variant of an earlier program called RADSPLIERE (Ref. 14) for 
computing the radiation from an internally-driven submerged spherical shell. 


Numerical Solution 

A NASTRAN finite element model of the spherical shell was prepared using 40 
axisymmetric conical shell elements spanning the 180 degrees between the two poles of the 
sphere. Due to the axisymmetry of the incident pressure loading, only the N ~ 0 harmonic 
was required. Since all structural points are in contact with both interior and exterior fluids, 
the resulting model therefore had 205 independent structural degrees of freedom (DOF) and 
41 fluid DOF for each of the two fluid domains. System matrices for the exterior fluid were 
also augmented by the addition of four constraint equations associated with interior Chief 
points to ensure uniqueness of the integral representation at the upper frequencies. The 
nondimensional frequency range 0<ka<5 was swept using a frequency increment of about ka 
= 0.05 with NASHUA and ka = 0.005 with the series solution. Since the series solution is 
converged, we treat it as an ’’exact" solution for this problem. 

The comparison between the computed and exact solutions is presented is Figs. 3 and 
4, which plot the frequency response of the nondimensional scattered pressure pr/(p Q a), 
where p is the far-field scattered pressure at distance r from the origin, p Q is the incident 
pressure, and a is the mean radius of the spherical shell. These two figures show very good 
agreement between the two scattering solutions in the backward (0 = 0) and forward ( 0 = 180 
degrees) directions. In fact, the computed and series solutions are virtually indistinguishable 
from each other. 
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Fig, 3, Forward Scattering from a Fluid-Filled Spherical Shell 



K A 


Fig. 4. Backward Scattering from a Fluid-Filled Spherical Shell 

DISCUSSION 

A very general computational capability has been described for predicting the sound 
pressure field radiated or scattered by arbitrary, submerged, fluid-filled, three-dimensional 
elastic structures subjected to lime-harmonic loads. The structure is modeled with 
NASTRAN (in all the generality that NASTRAN allows) in combination with boundary 


162 





element models of both interior and exterior fluid domains. Sufficient automation is 
provided so that, for many structures of practical interest, an existing structural model can 
be adapted for NASHUA acoustic analysis within a few hours. 

One of the many benefits of having NASHUA linked with NASTRAN is the ability to 
integrate the acoustic analysis of a structure with other dynamic analyses. Thus the same 
finite element model can be used for modal analysis, frequency response analysis, linear 
shock analysis, and underwater acoustic analysis. In addition, many of the pre- and 
postprocessors developed for use with NASTRAN become available for NASHUA as well. 
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